{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example-3-MC-GPA-BS-CVA\n",
    "# Author: Matthew Dixon\n",
    "# License: MIT\n",
    "# Email: matthew.dixon@iit.edu\n",
    "# Notes: tested on Mac OS X running Python 3.6.9 with the following packages:\n",
    "# scikit-learn=0.22.1, numpy=1.18.1, matplotlib=3.1.3\n",
    "# Citation: Please cite the following reference if this notebook is used for research purposes:\n",
    "# M.F. Dixon, I. Halperin and P. Bilokon, Machine Learning in Finance: From Theory to Practice, Springer Graduate textbook Series, 2020. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview\n",
    "The purpose of this notebook is to demonstrate the use of a Gaussian Process Regression model (GP) in CVA modeling of a derivative portfolio. In this notebook, European option prices are generated from the Black-Scholes model.  The notebook illustrate the CVA estimation for a portfolio which is short a put and long two calls.  The GP model is combined with an Euler time-stepper to generate a GP MtM cube. The expected positive exposure of the portfolio is compared using exact derivative prices and GP derivative prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn import gaussian_process\n",
    "from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, RBF\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')\n",
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "import sys\n",
    "sys.path.append('../')\n",
    "from BlackScholes import bsformula\n",
    "from GBMSim import gbm\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train GPs to a static representation of the option price surface at time snapshots of maturities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainGPs(x_train, f, timegrid):\n",
    "\n",
    "    gps = []\n",
    "    i = 0\n",
    "    for time in timegrid:\n",
    "        \n",
    "        y_train = []\n",
    "        \n",
    "        for idx in range(len(x_train)):\n",
    "             y_train.append(f(x_train[idx], time))\n",
    "        y_train = np.array(y_train)\n",
    "\n",
    "        sk_kernel = RBF(length_scale=1.0, length_scale_bounds=(0.01, 100.0))  #100000.0\n",
    "        gp = gaussian_process.GaussianProcessRegressor(kernel=sk_kernel, n_restarts_optimizer=20)\n",
    "        gp.fit(x_train,y_train)\n",
    "        gps.append(gp)\n",
    "        i += 1\n",
    "    return gps   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict the option prices over a timegrid of different maturities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predictGPs(x_test, f, gps, timegrid):\n",
    "    \n",
    "    preds = []\n",
    "    stds = []\n",
    "    y_tests = []\n",
    "    i = 0\n",
    "    for time in timegrid:\n",
    "        \n",
    "        y_test_= []\n",
    "        for idx in range(len(x_test)):\n",
    "            y_test_.append(f(x_test[idx], time))\n",
    "            \n",
    "        y_test = np.array(y_test_)\n",
    "        y_tests.append(y_test)\n",
    "        \n",
    "        # Find optimal model hyperparameters\n",
    "        # Set into eval mode\n",
    "        pred, std = gps[i].predict(x_test,return_std=True)\n",
    "        preds.append(pred)\n",
    "        stds.append(std)\n",
    "        i+=1\n",
    "    return y_tests, preds, stds "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Black-Scholes derivative portfolio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify the portfolio\n",
    "KC = 110    # Call strike\n",
    "KP = 90     # Put strike\n",
    "lb = 0.01\n",
    "ub = 400\n",
    "r  = 0\n",
    "sigma = 0.3\n",
    "nt = 11\n",
    "T = 2.0\n",
    "S0 = 100\n",
    "timegrid = np.array(np.linspace(0.0 ,T, nt), dtype='float32').reshape(nt, 1)\n",
    "\n",
    "\n",
    "portfolio = {}\n",
    "portfolio['call'] = {}\n",
    "portfolio['put'] = {}\n",
    "\n",
    "portfolio['call']['price']= lambda x,y: bsformula(1, lb+(ub-lb)*x, KC, r, y, sigma, 0)[0]\n",
    "portfolio['put']['price']= lambda x,y: bsformula(-1, lb+(ub-lb)*x, KP, r, y, sigma, 0)[0]\n",
    "portfolio['call']['weight']=2.0\n",
    "portfolio['put']['weight']=-1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare training and test data\n",
    "training_number= 100\n",
    "testing_number = 50\n",
    "x_train = np.array(np.linspace(0.0,1.0, training_number), dtype='float32').reshape(training_number, 1)\n",
    "x_test = np.array(np.linspace(0.0,1.0, testing_number), dtype='float32').reshape(testing_number, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train and predict the portfolio of derivative prices over risk horizon.\n",
    "for key in portfolio.keys():\n",
    "    portfolio[key]['GPs'] = trainGPs(x_train, portfolio[key]['price'], timegrid)\n",
    "    portfolio[key]['y_tests'], portfolio[key]['preds'], portfolio[key]['sigmas'] = predictGPs(x_test, portfolio[key]['price'], portfolio[key]['GPs'], timegrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CVA simulation using the credit default model described in the above reference.\n",
    "def CVA_simulation(sim_params, model_params, def_model):\n",
    "    \n",
    "    n_sim_dt = sim_params['n_sim_dt'] # number of Euler stpes\n",
    "    M        = sim_params['M']        # number of paths\n",
    "    nt       = sim_params['nt']       # number of exposure dates\n",
    "    timegrid = sim_params['timegrid'] # time grid of exposure dates\n",
    "    \n",
    "    r        = model_params['r']\n",
    "    sigma    = model_params['sigma']\n",
    "    T        = model_params['T'] \n",
    "    t0       = model_params['t0'] \n",
    "    S0       = model_params['S0'] \n",
    "    \n",
    "    gamma_0  = def_model['gamma_0']\n",
    "    gamma_1  = def_model['gamma_1']\n",
    "    \n",
    "    \n",
    "    stride = n_sim_dt/(nt-1)\n",
    "    idx = np.arange(0,n_sim_dt+1,stride, dtype=int)\n",
    "    \n",
    "    pi = {}\n",
    "    pi['tilde'] = np.array([0.0]*(nt-1)*M, dtype='float32').reshape((nt-1), M)     # GP portfolio value\n",
    "    pi['exact'] = np.array([0.0]*(nt-1)*M, dtype='float32').reshape((nt-1), M)     # BS portfolio value\n",
    "    pi['tilde_var'] = np.array([0.0]*(nt-1)*M, dtype='float32').reshape((nt-1), M) # GP portfolio variance\n",
    "    gamma = np.array([0.0]*(nt-1)*M, dtype='float32').reshape((nt-1), M)           # hazard rates\n",
    "    dPD = np.array([0.0]*(nt-1)*M, dtype='float32').reshape((nt-1), M)             # default probabilities\n",
    "    \n",
    "    #simulate underlying GBM dynamics using Euler\n",
    "    S = gbm(S0, r, sigma, T-t0, n_sim_dt, M)\n",
    "    \n",
    "    if (def_model['calibrate']):\n",
    "        x  = np.exp(S0/S)**gamma_1\n",
    "        # default probability (assumed to be estimated from credit spread)\n",
    "        dt = timegrid[1]-timegrid[0] \n",
    "        f  = lambda y: np.abs(np.mean(np.prod(x**(-y*dt), axis=0)) - def_model['p']) \n",
    "        res = sp.optimize.basinhopping(f, 0.1, niter=10)\n",
    "        i = 1   \n",
    "        while (abs(res.fun) >1e-3):   \n",
    "          res = sp.optimize.basinhopping(f, 0.1, niter=100*i)\n",
    "          i *= 2\n",
    "        gamma_0= res.x[0]\n",
    "        print(\"calibration:\", gamma_0, gamma_1, f(gamma_0), res.fun)  \n",
    "    \n",
    "    \n",
    "    for m in range(M):  \n",
    "      i = 1 \n",
    "      exp_factor=1\n",
    "        \n",
    "      for time in timegrid[1:]:\n",
    "        dt = timegrid[i]-timegrid[i-1] \n",
    "        \n",
    "        S_= S[idx[i],m] # simulated S\n",
    "        # avoid simulated S breaching boundaries of domain\n",
    "        if (S_<lb):\n",
    "            mins=S_\n",
    "            S_=lb\n",
    "        if (S_>ub):\n",
    "            S_=ub\n",
    "            maxs=S_\n",
    "    \n",
    "        pred_= 0\n",
    "        v_ = 0\n",
    "        var_ =0 \n",
    "    \n",
    "        for key in portfolio.keys():\n",
    "           pred, std = portfolio[key]['GPs'][i].predict(np.array([(S_-lb)/(ub-lb)]).reshape(1,-1),return_std=True) \n",
    "           pred_ += portfolio[key]['weight']*pred\n",
    "           var_ += (portfolio[key]['weight']*std)**2 \n",
    "       \n",
    "           if key=='call':\n",
    "              v_ += portfolio[key]['weight']*bsformula(1, S_, KC, r, time, sigma, 0)[0]\n",
    "           else:\n",
    "              v_ += portfolio[key]['weight']*bsformula(-1, S_, KP, r, time, sigma, 0)[0]\n",
    "        pi['tilde'][i-1,m] = np.maximum(pred_,0)\n",
    "        pi['exact'][i-1,m] = np.maximum(v_,0)\n",
    "        pi['tilde_var'][i-1,m] =var_ \n",
    "          \n",
    "        # default intensity model\n",
    "        gamma[i-1,m] = gamma_0*(S0/S_)**gamma_1    \n",
    "        \n",
    "        # compute default probabilities  \n",
    "        exp_factor*=np.exp(-dt*gamma[i-1,m])    \n",
    "        dPD[i-1,m]= gamma[i-1,m]*exp_factor\n",
    "        \n",
    "        i += 1\n",
    "    # compute CVA\n",
    "    i = 0\n",
    "    CVA ={}\n",
    "    CVA['tilde'] = 0\n",
    "    CVA['exact'] = 0\n",
    "    CVA['tilde_up'] = 0\n",
    "    CVA['tilde_down'] = 0\n",
    "   \n",
    "    for time in timegrid[1:]:\n",
    "        dt = timegrid[i+1]-timegrid[i]\n",
    "        mu_tilde = np.mean(dPD[i,:]*pi['tilde'][i,:])*np.exp(-r*(time-t0))*dt\n",
    "        CVA['tilde'] += mu_tilde\n",
    "        std_err_MC = np.std(dPD[i,:]*pi['tilde'][i,:])*np.exp(-r*(time-t0))*dt/np.sqrt(M)\n",
    "        CVA['tilde_up'] += mu_tilde + 2.0*std_err_MC\n",
    "        CVA['tilde_down'] += mu_tilde - 2.0*std_err_MC\n",
    "        CVA['exact'] += np.mean(dPD[i,:]*pi['exact'][i,:])*np.exp(-r*(time-t0))*dt\n",
    "        i+=1\n",
    "    CVA['tilde'] *= (1 - def_model['recovery'])\n",
    "    CVA['tilde_up'] *= (1 - def_model['recovery'])\n",
    "    CVA['tilde_down'] *= (1 - def_model['recovery'])\n",
    "    CVA['exact'] *= (1 - def_model['recovery'])\n",
    "        \n",
    "    return(CVA, pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_params = {}\n",
    "model_params = {}\n",
    "def_model = {}\n",
    "  \n",
    "model_params['r'] = r        # risk-free rate\n",
    "model_params['sigma'] = sigma  # implied volatility\n",
    "model_params['T'] = T        # Time to maturity \n",
    "model_params['t0'] = 0\n",
    "model_params['S0'] = S0     # Underlying spot\n",
    "     \n",
    "# parameters of the default intensity model    \n",
    "def_model['gamma_0'] = 0.01\n",
    "def_model['gamma_1'] = 1.2\n",
    "def_model['calibrate'] = False    \n",
    "def_model['recovery'] = 0.4\n",
    "\n",
    "sim_params['n_sim_dt'] = 100 # number of Euler stpes\n",
    "sim_params['M']  = 1000      # number of paths\n",
    "sim_params['nt'] = nt        # number of exposure dates\n",
    "\n",
    "sim_params['timegrid']= timegrid # time grid of exposure dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train and predict over portfolio, then perform the CVA simulation to give the GP MTM cube\n",
    "for key in portfolio.keys():\n",
    "    portfolio[key]['GPs'] = trainGPs(x_train, portfolio[key]['price'], timegrid)\n",
    "    portfolio[key]['y_tests'], portfolio[key]['preds'], portfolio[key]['sigmas'] = predictGPs(x_test, portfolio[key]['price'], portfolio[key]['GPs'], timegrid)\n",
    "\n",
    "CVA_0, pi = CVA_simulation(sim_params, model_params, def_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tilde': array([0.1689328], dtype=float32),\n",
       " 'exact': array([0.16893284], dtype=float32),\n",
       " 'tilde_up': array([0.1833201], dtype=float32),\n",
       " 'tilde_down': array([0.15454549], dtype=float32)}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compare the exact CVA, i.e. using BS prices, with the GP portfolio and provide the GP error bounds \n",
    "CVA_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x119389828>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnkAAAF2CAYAAAASt47pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdfXxT9d0//tfJOWl636ZpoZaWKSg6BlO0iIBy14IKXsoURZ06dZfoUEG8kO9wc+ocdxcyGHLjDRWnXo8NHY4xnYroBig6YYCIuAcggkDvkiZt0zRNcnLO7w9+OSM0aU/bnCZtXs/HgwfNyTnJOyd377w/d4KqqiqIiIiIqFcxxTsAIiIiIoo9JnlEREREvRCTPCIiIqJeiEkeERERUS/EJI+IiIioF2KSR0RERNQLMckjIiIi6oWkeAeQiCorK+MdQlT5+flwOBzxDiPh8LxEx3MTGc9LdDw3kfG8RMdzE1l3nJeioqKo17GSR0RERNQLMckjIiIi6oWY5BERERH1QuyTR0REMaeqKlpaWqAoCgRBiHc4MVFTUwOfzxfvMBISz01ksTovqqrCZDIhNTW1Q+8nJnlERBRzLS0tMJvNkKTe8zUjSRJEUYx3GAmJ5yayWJ4XWZbR0tKCtLQ03cewuZaIiGJOUZReleARxZskSVAUpUPHMMkjIqKY6y1NtESJpKPvKyZ5RETUK5WUlGDixInav1WrVhl2X1988QWeeOKJDh0zYsQIlJWVafF19PjOWrlyZYeP2bBhA37xi19E3D506FBMmDAB48aNw//93/91+Larq6tx3333AQAOHDiADz/8ULtuy5YtMXneQnGe+Xo4dOhQl2830bGWTkREvVJqaio++OCDNvcJBoNhfaZkWdbVzHz2fhdffDEuvvjiDsf45ptvIi8vr8PHdcVzzz2HWbNmxez2rr/+eixZsgTV1dUYP348Jk2ahIKCAt3HFxYW4qWXXgIAfPXVV9i/fz/KysoAAJMmTcKkSZNiFueCBQticluxpPc11xms5BERUVIZMWIEli9fjqlTp+Ltt9/GtGnTsGjRItx0001Yt24dTp48iVtuuQXl5eW45ZZbcOrUKQDArFmz8NRTT2HatGmtkoWdO3firrvuAgAsW7YMjz76KKZNm4aRI0eioqJCd2yyLGPy5MnYuXMnAGDRokVYvHixFveCBQswZcoUTJkyBd9++y0AoK6uDvfddx8mT56MyZMnY9euXQAAj8eDOXPmoKysDOXl5XjnnXewcOFCtLS0YOLEiXjooYcAABs3bsSUKVMwceJEzJs3D8FgEMDp6teVV16Jm266Cbt372439vz8fHzve9/DyZMn4XK5cO+996K8vBzXXXcdDh48CAD49NNPtUrapEmT0NTUhBMnTmDChAnw+/149tlnsXnzZkycOBF/+ctftApiY2MjRowYofVJ83q9KC0tRSAQwLFjx/DjH/8Y11xzDX70ox/hyJEjus/3u+++i+nTp0NVVdTU1ODKK69EbW0tNmzYgHvuuQc//vGPcdVVV+G3v/2tdswLL7yACRMmYMKECVpy2tzcjDvvvBPl5eWYMGEC/vKXvwAASktL4XQ6AZyu9k6bNg3A6dfIvHnzcNttt2H27NkIBoN45plnMHnyZJSXl+O1117T/RjawkoeEREZ6le/+pX2JR8rgwcPxq9//es29wklMyEPPfQQbrjhBgCAxWLBpk2bAACvvfYaGhsbsXHjRgDAT37yE0ybNg233HIL/vjHP+KJJ57Ayy+/DAA4evQoNmzY0O6IySNHjuDNN9+Ex+PBVVddhbvuugtms7nVfjfffDNMJpP294wZM7B8+XLMmDEDzzzzDP7+97/j7bff1vbPzMzEO++8gzfffBNPPvkkXn31VfzqV7/Cfffdh8svvxynTp3C7bffjm3btmHFihXIysrSmj/r6+sxZcoUrF+/XqtwHj58GJs3b8amTZtgNpsxf/58vPXWWxgzZgyeffZZvPfee8jKysLNN9+MIUOGtPmYjx8/ju+++w7nnnsuli1bhiFDhuDll1/Gxx9/jNmzZ+ODDz7A888/j4ULF2L48OHweDywWCza8SkpKZg7dy7279+vJdEbNmwAAGRnZ2Pw4MH49NNPMXr0aGzZsgXjxo2D2WzGvHnzsHjxYgwYMAB79uzB/Pnz8eabb7aKb/Pmzfj888/DLl977bX429/+hldeeQV///vfMXfuXPTp0wcAsG/fPnz44YdIS0vDlClTUFZWBkEQ8MYbb+Dtt9+Gqqq47rrrMHLkSBw/fhyFhYVactbY2NjmuQKA/fv3489//jPS0tLw+uuvIysrC3/729/g8/kwdepUjB07Fv3792/3dtqSMEnevn37sH79eiiKgrKyMkydOjXselVVsX79euzduxcWiwUzZ87EgAEDAABr1qzBnj17kJOTg2XLlmnHvPHGG/jwww+RnZ0NALjttttw6aWXdt+DIiKiuGmrufb666+Pevlf//oX1q1bBwC46aab8Jvf/Ea77rrrrtM1JUZZWRksFgssFgvy8/Nht9sjrjEaqbn2wgsvxE033YS7774bmzdvRkpKinZd6Ltx6tSpeOqppwAAO3bsCOtf1tTUhKamJuzYsQNr1qzRtufm5ra6/48//hhffvklJk+eDOB0Ypyfn4+9e/di5MiRsNls2vk5evRoxMe6efNm7Nq1CykpKViyZAmsVis+//xzrcp15ZVXwuVyobGxEcOHD8fTTz+NH/3oR7j22mvbXHf1bNdffz02b96M0aNHY/PmzfjJT34Cj8eDf/3rX7j//vu1/fx+f9TjIzXXPvPMMygrK8Oll14alntcddVV2nNz7bXX4vPPP4cgCLjmmmuQnp6ubf/nP/+JcePG4ZlnnsGCBQtQXl6OESNGtPt4Jk2apE2Hsm3bNnz99dd45513AAButxvffvtt70jyFEVBRUUFfvnLX8Jms2H+/PkoLS1FcXGxts/evXtRXV2NlStX4vDhw1i3bh0WLlwIABg3bhyuueYarF69utVtT5kypdWbmYiIuk97Fbd4CH1JR7t8pjNHNLa135nOrFCJoqg1ger173//G9nZ2bDb7VFjCf2tKAo2b97cav40VVXbHY2pqipuvvlmzJ8/P2z7e++9p3skZ6hPnizLYbd7NkEQ8NBDD6GsrAwfffQR/uu//gt/+MMfkJqaqut+Jk2ahEWLFsHlcmH//v0YPXo0mpubkZ2d3W7fy7ZUV1dDEATY7XYoiqJVVs9+/IIgRHxcADBw4EC8++67+Oijj7Bo0SKMHTsWc+bMgSiKWhPz2ZMin/1a+s1vfoNx48Z1+nFEkhB98o4cOYLCwkL07dsXkiRh1KhRWp+CkN27d2PMmDEQBAGDBg2Cx+OBy+UCcLpsn5mZGY/QiYiolyktLdX6VL311lu4/PLLu/X+//a3v8HlcmHjxo144okn0NDQoF23efNm7f/LLrsMADB27Fi88sor2j4HDhzQtq9fv17bXl9fDwAwm80IBAIATlfZ3n77bTgcDgCAy+XCyZMnMWzYMHz66adwOp0IBAJhTcZ6XHHFFXjrrbcAnO6vmJeXh6ysLBw7dgzf//738eCDD2Lo0KH497//Db/fD1VV4ff7kZqaisbGRgSDwVZzwmVkZOCSSy7Br371K5SXl0MURWRlZaGkpAR//etfAZxOLr/66ivdccqyjEcffRSrV6/GBRdcgBdffFG7bseOHXC5XPB6vXj//fcxfPhwXHHFFXj//ffh9XrR3NyM9957DyNGjEB1dTXS0tJw00034YEHHsCXX34J4PQI7/379wOAVqWLZOzYsXj11Ve15+Wbb75Bc3Oz7scRTUJU8pxOp1YSBgCbzYbDhw+32ic/Pz9sH6fTCavV2uZtv//++9i+fTsGDBiAu+66K2IyuHXrVmzduhUAsHjx4rD7STSSJCV0fPHC8xIdz01kPC/RxeLc1NTUxH0y5JaWlrCRmePHj8cTTzwBQRAgiqIW39mXFy5ciEceeQTPP/88bDYbfve732nXnbnfmURRhCAIkCQJJpMJJpMpbL9IxwmCgJtvvllr/h08eDCefvppLFq0CH/605/Qr18//PSnP8VTTz2F5557DoIgQJZlXHfddVBVFc8//zwkScLChQvx85//HOXl5QgGg7jiiiuwdOlS/M///A9+/vOfY8KECRBFEXPnzsWUKVNw5513YuLEiRg6dCjWrl2L+fPn4/bbb4eiKDCbzVi0aBFKS0vx2GOP4YYbbkDfvn3xwx/+MOIE16IoapWvM6+bN28eZs+ejfLycqSlpeG5556DJEmoqKjAJ598AlEUMXDgQJSVlaG2thYAYDKZMHr0aDz//PO4+uqrMXPmTC3ZCwaDEAQB119/Pe6//3689dZb2v2tXbsW/+///T+sXLkSsixj6tSprUY6i6KIv/71r2EFpCVLlmD79u244oorMHr0aFx88cW4+uqrMWnSJIiiiBEjRuCRRx7Bt99+ixtvvFFLqm+99VZMmTIFAPDjH/8Yl1xyCf7+97/j6aefhslkgtlsxpIlSyBJEubOnYs5c+Zg1apVuPTSS6O+Ru666y6cOnUK11xzDVRVhc1mw+9///tW5zvU/K+XoEarPXajTz/9FF988QUeeOABAMD27dtx5MgR3Hvvvdo+ixYtwo9+9CNcdNFFAE6X/++44w6tX15tbS2WLFkS1ievvr5e64+3YcMGuFwuzJw5s914KisrY/bYYi0/P1/7xUX/wfMSHc9NZDwv0cXi3DQ3N+tu2uwpJEkKa5LsbiNGjMC7777b7VOu6NGRcxOq2nVWKG0JNacKggBBEGAymbS/uzoZ94YNG8IGgHRWrF8zkd5XbfVrTIjmWpvNhrq6Ou1yXV1dqwqdzWYL+9CJtM/ZcnNztWy5rKwM33zzTWwDJyKiuEmAGgV1QjAY7NJzd3YSp6oqFEVBIBCA3++Hz+eDz+eD3++H3++HLMvafSbbayYhkryBAweiqqoKtbW1kGUZO3fuRGlpadg+paWl2L59O1RVxaFDh5Cent5ukhfqswcAn3/+OUpKSgyJn4iIupeiKPD5fEn1pf3Pf/4zIat4HaGqKmRZNmTZu7OreKGkLhgMIhAIaMlfKAEMBAKQZRmKorR6HU2fPj0hJ07uqITokyeKIu69914sWLAAiqJg/PjxKCkpwZYtWwCcHlEzbNgw7NmzB7NmzUJKSkpYs+uKFStw8OBBuN1uPPDAA7jlllswYcIEvP766zh27BgEQUBBQQFmzJgRr4dIREQxFEoUZFmOOP8cJaZAIBCXdY3Pvs9QAnhmcndmghjr5t94SYg+eYmGffJ6Hp6X6HhuIuN5iS4W58bj8SAjIyNGEYVTFAV+v1+b0iIlJUXr/G+kePfJS2R6zk2ootZTEqZICWCk/n9tifVrJtL7qq0+eQlRySMiot7FZDIZtibnmc19giAgEAggJSWlxyQPycjIZlqjRKv+AdASt+6s/smy3OEfM0zyiIgo5lJTU9HS0gKfzxfTL7tgMAiPxxP2ZaeqKiRJajUZcKxZLJZWE9rSae2dm5aWlh5VxeusM+f2M5lMSE9Pj0mVWVVVmEwm3RNHhzDJIyKimBMEwZCk69SpUxGThUAggMzMzLCVJmKNTfzRtXVu/H4/nE5nUvadTE1NjetiDQkxupaIiKg9Pp8PXq83YjVIkiRUVVUl1WjbnqK6ujruE2MnKyZ5RETUI9jt9qjJgiAICAaDYXOuUvw1NDQkRTNtomKSR0RECa+lpSVqFS9EkiQ0NDSw31yCCAaDcDgcrOLFEZM8IiJKeHa7XVefLlEUUV1dzWbbBFBbW9stU9tQdDz7RESU0Lxer+5RuqFmW6fT2Q2RUTTNzc2tRkFT9+PZJyKihNbRJj9RFOFyueD3+w2MiqJRVRW1tbVJOZo20TDJIyKihNXc3IyWlpYOd9yXJAnV1dUGRUVtqaurQzAYjHcYBCZ5RESUwBwOR6cqQqGVMFwulwFRUTR+vx/19fUcbJEgmOQREVFCam5u1tao7QxJkuB0OrnebDfinHiJhUkeERElJL0jatsSGm1LxuOceImHSR4RESUcj8eDQCDQ5dsRBAE+nw8NDQ0xiIqiCU1EzSpeYmGSR0RECaezffEikSQJDoeDgwEMVFlZyQpeAmKSR0RECcXtdse8H53JZGKzrUG8Xi8aGxs5J14C4jNCREQJQ1VVQ5r9TCYTvF4v3G53TG832amqipqaGqSkpMQ7FIqASR4RESWMpqYmw0bDms1m2O12KIpiyO0nI6fTyWbwBMYkj4iIEkKoimfkSgmCIKCmpsaw208mgUCAc+IlOCZ5RESUEBobGw2vCplMJng8Hng8HkPvJxnU1NRAFMV4h0FtYJJHRERxp6oqXC5Xt1SFzGYzamtr2WzbBW63Gz6fjyNqExyTPCIiirvuqOKdSVVV1NbWdtv99SaKosBut7OZtgdgkkdERHGlqiqcTme3Jg2iKKKpqQler7fb7rO3qK2tZQWvh2CSR0REcdXQ0ABVVbv9fs1mM2pqauJy3z2V1+tFU1MT58TrIfgsERFR3IT64sWrA3+o6ZHaF2riZjNtz8Ekj4iI4qa+vj6uAyBEUURjYyN8Pl/cYugpQnPisam252CSR0REcdGdI2rbIkkSqqqq2GzbBlmW41pxpc5hkkdERHHhcrkSIrESBAHBYBB1dXXxDiVh1dTUxD0Zp45jkkdERN1OUZSEWi1BkiTU19ez2TYCt9sNr9fLZtoeiEkeERF1O5fLFe8QWpEkiaNtzxIamGLkUnNkHCZ5RETUrRRFQUNDQ8L17xIEQet7Rqc5HI54h0BdwCSPiIi6ldPpjHcIUYmiCJfLhUAgEO9Q4s7n86GxsTHhknHSj0keERF1m0St4p1JFEVUVVXFO4y4UlUV1dXVCdNnkjqHSR4REXWburq6hO/ALwgCAoFAUjfb1tfXc068XoBJHhERdQtFUXpM858kSXA6nZBlOd6hdDtZluF0OnvE80Rt61CSpyhKUv+yISKiznM4HD1qzVNRFFFdXR3vMLpdTU1Nj3qeKDpdje0ejwfr1q3DZ599BkmS8Nprr2H37t04cuQIbr31VqNjJCKiHi4YDMLtdveoPl6CIMDn86GhoQE5OTnxDqdbeDweeL1eTpnSS+hK1V966SWkp6djzZo12ht00KBB2Llzp6HBERFR71BXV9cjq0OSJMHhcCAYDMY7FMMpioLa2lomeL2Irnfcl19+iXvuuQdWq1Xblp2djYaGBsMCIyKi3iEYDKKxsbFHJnkAYDKZUFNTE+8wDOdwODgRdC+jq26enp4Ot9sdluQ5HI6wy121b98+rF+/HoqioKysDFOnTg27XlVVrF+/Hnv37oXFYsHMmTMxYMAAAMCaNWuwZ88e5OTkYNmyZdoxTU1NWL58Oex2OwoKCjBnzhxkZmbGLGYiImqfw+Ho0Z34TSYTmpube3VhIzQnHqt4vYuun1VlZWVYtmwZDhw4AFVVcejQIaxevRoTJ06MSRCKoqCiogKPP/44li9fjk8++QQnT54M22fv3r2orq7GypUrMWPGDKxbt067bty4cXj88cdb3e6mTZswdOhQrFy5EkOHDsWmTZtiEi8REekjyzLcbnePreKFmM1mVFZWQlGUeIcSc6qqoqampkf1lyR9dL3rbrjhBowcORIVFRUIBoNYu3YtSktLMXny5JgEceTIERQWFqJv376QJAmjRo3Crl27wvbZvXs3xowZA0EQMGjQIHg8Hm2k7+DBgyNW6Hbt2oWxY8cCAMaOHdvqNqnrgsEgvF4v6uvrUVNTg5MnT+L48eM4evQoF/omoh5fxTtTb222bWhoQCAQ4Jx4vVC7abuiKPjHP/6BSZMmYcqUKYYE4XQ6YbPZtMs2mw2HDx9utU9+fn7YPk6ns80m44aGBu16q9WKxsbGGEeeHBRFgSzL8Pl88Hq9kGVZ+xf6VWsymWAymbQPCZPJhJMnT6J///4s/xMlKVmW0dTU1Gs+A0wmEzweDzweDzIyMuIdTkwEg0HU1dWxitdLtfusmkwmvPrqq5gwYYJhQUTq6Hn2Lwo9+3TW1q1bsXXrVgDA4sWLw5LJRCNJkiHxqaoKWZbh9/vh8Xjg8/kQCATCkjmTyQRJknR/YKuqioaGBpx//vmGf4AYdV56A56byHheoovVufnuu+9gtVp7fFNtiCiKsNls8Hq9KCkp6RWP6/jx48jNze3yYxFFEdnZ2TGKqvcQRTGunzO6vnkvu+wy7N69G6WlpYYEYbPZUFdXp12uq6trVaGz2WxwOBxt7nO2nJwcuFwuWK1WuFyuqC/A8vJylJeXa5fPvJ9Ek5+f36X4gsEg/H4/fD6flsgFg0HIsqwl0qIoxuzDS1VVfPHFF4Z/IHb1vPRmPDeR8bxEF4tzI8syKisre00VDzg9q0RjYyOCwSAOHjyIwsLCeIfUJR6PB1VVVTF5jkLnhsJZrVbDP2eKioqiXqcryQsEAvjtb3+LQYMGwWazhVXQHnrooS4HOHDgQFRVVaG2thZ5eXnYuXMnZs2aFbZPaWkp3nvvPYwePRqHDx9Genp6u0leaWkptm3bhqlTp2Lbtm0YPnx4l2PtCULNq36/H16vV0vkAoFA1OZVoyptgiBAURScPHkSJSUl7PNBlCRqa2t7bROgKIpoamqC1+tFWlpavMPpFFVVYbfbe1USTq3pegeWlJSgpKTEsCBEUcS9996LBQsWQFEUjB8/HiUlJdiyZQsAYNKkSRg2bBj27NmDWbNmISUlBTNnztSOX7FiBQ4ePAi3240HHngAt9xyCyZMmICpU6di+fLl+Oijj5Cfn49HH33UsMfQ3VRV1apyoUTu7KqcyWSCKIpaYiWKYlw6QJtMJu1XfVFRERM9ol7O7/ejubm5VycQZrMZ1dXVOPfcc3vkZ5rD4YCiKL1mUAxFJqic+bCVysrKeIegCSVtXq8XPp8P6enp2qLZRjSvGkmWZaSnp+Occ86J+W2z6S06npvIeF6i6+q5OXXqVK8crXl2k2QwGERmZib69OkTx6g6zu/347vvvotpEs7m2sisVqvhfRW73Fx74MCBqNcNGTKk4xFRGFVVEQgEojavqqqqJXIpKSkAjGteNZIkSWhubkZtbW2P+1AkIn1Co/B7cxUvRBRFNDY2IicnBxaLJd7h6FZdXd0jv0Oo43Q9y2vXrg273NjYCFmWYbPZsGrVKkMC681CE08mavOqkSRJgtvt1kapEVHv4nA4kiqBkCQJ1dXV6N+/f4+oXDY0NMDv9ydFEk46k7zVq1eHXVYUBRs3buyxHU7jTVVVbe4oQRCS7s0mSRLq6+shSRJycnLiHQ4RxUgyVfFCBEGALMut5ntNRMFgEA6HI6men2TXqY5cJpMJN954I/7yl7/EOh5KEpIkwW63w+12xzsUIooRu92eVFW8EEmS4HK54Pf74x1Km2pra3tE/22KnU4/2/v37+eLhbrEbDajpqYGzc3N8Q6FiLqopaUFXq+3RzRZGiHUbJuoYxmbm5vh8Xj4vZ1kdP3k+tnPfhZ22e/3w+/347//+78NCYqSh9lsRlVVFYqLi3tUx2UiCpfsc66Fmm1dLhfy8vLiHU6YUD/wZH5+kpWuJO/hhx8Ou2yxWHDOOecgPT3dkKAouYiiiFOnTqG4uFgbPUxEPUdoiqdkTyJEUYTL5UJWVlZCnYu6ujptaUpKLrqSvMGDB4dd9vv9fLFQzAiCAJPJhJMnT6J///5J2aeHqCdLthG1bRFFEVVVVejfv3+8QwFw+vu6vr4+oZJO6j66MrVXX30VR44cAQDs2bMH99xzD+6++27s3r3b0OAoeZyZ6IWWXiOixNfc3IyWlpak7Yt3NkEQEAgEUF9fH+9QAHBOvGSnK8n7+OOPtWXN/vSnP+Hhhx/GvHnz8Ic//MHQ4Ci5CIIAVVVx4sQJJnpEPQSn5GhNkiTU1dVBluW4xtHQ0NArVx4h/XQleT6fDxaLBW63GzU1Nbjiiivwwx/+kEsCUcyZTCYEg0FUVlYm7Cg1IjqtubkZfr+fSUQEoiiipqYmbvevKArq6upYxUtyup79oqIi7NixA9XV1fjhD38I4PSqF+wkT0YQRRE+nw/V1dUoLCzkFwhRgkr2EbVtEQQBLS0taGhoiMuk7zU1NfzsJH2VvJ/+9Kd4//33ceDAAUyfPh0A8MUXX2gJH1GsSZIEr9eL2traeIdCRBF4PB4EAoF4h5HQJEmCw+FAMBjs1vv1er2cE48A6KzknX/++fjNb34Ttu2qq67CVVddZUhQRMDpil5TUxMkSUr45YKIkg374uljMplQU1ODoqKibrk/zolHZ9LdWH/gwAFs374dLpcLVqsVY8aMwZAhQ4yMjUhb59ZkMsFqtcY7HCIC4Ha7Icsy+3vpYDKZ0NzcjKamJmRmZhp+f06nE8FgkM8NAdDZXPvhhx9ixYoVyM3NxeWXXw6r1Yrf/e532Lp1q9HxEWkj1RobG+MdClHSU1WVHfo7yGw2o7a21vBZA2RZRn19PZ8b0uh6JWzevBm//OUvce6552rbRo0ahWXLlqG8vNyo2Ig0oQ9JURSRkZER73CIklZTUxNkWWZzYAcJgoCamhqcc845ht1HdXU1RFE07Pap59FVyXO73SguLg7bVlRUhKamJkOCIooktM6t1+uNdyhESSlUxWOC13EmkwkejwfNzc2G3L7b7YbP5+OIWgqjK8m76KKL8Oqrr8Ln8wEAWlpa8Nprr2HQoEGGBkd0NrPZjMrKSvj9/niHQpR0Ghsbu32kaG9iNptRU1MT82ZbRVFgt9vZTEut6HpF3HfffVixYgXuvvtuZGZmoqmpCYMGDcLs2bONjo+oFVEUuc4tUTdTVRUul4vvuS5SVRV2ux19+/aN2W3W1taygkcR6Xq3Wq1WPP3006irq9NG13JKC4qX0Dq3J06cQP/+/dkHhagbhKp4TPK6RhRFuN1uZGdnIy0trcu319LSgqamJjahU0S6Z0r0eDw4ePCg9s/j8RgZF1GbQr9aT548yXVuiQymqiqcTicTvBgxm82orq7u8tKNoTnx+LxQNLqSvAMHDuDBBx/Eu+++iyNHjs+cQacAACAASURBVOC9997Dgw8+iC+//NLo+IiiMplMUBQFp06d4jq3RAZqaGjgeyzGVFXt8vrvLpcLwWCQTbUUla70v6KiAjNmzMCoUaO0bZ9++ikqKiqwYsUKw4Ijao/JZEIgEMDx48eRkZHBDzuiGAv1xWO3iNgSRRENDQ3Izs6GxWLp8PGyLLOPJLVLVyXP5XLhiiuuCNt2+eWXo76+3pCgiDpCFEV4vV7U1NTEOxSiXqe+vp5dIgwiSVKnm21ramqYeFO7dCV5Y8aMwXvvvRe2bcuWLRgzZowhQRF1lCRJaG5uht1uj3coRL0GR9QaSxAEyLIMp9PZoePcbje8Xi9bLqhdut653377LT744ANs3rwZeXl5cDqdaGhowAUXXIAnn3xS2+/pp582LFCi9oSaPyRJ4jq3RDHgcrnYF89gkiTB5XIhKysLKSkp7e6vKAocDgdH05IuupK8srIylJWVGR0LUZeZzWY4nU6Ioojs7Ox4h0PUYymKwnVQu4kkSaiqqkL//v3brc45HA4m3qSbrnfvuHHjDA6DKHYkSUJtbS1MJhMyMzPjHQ5Rj+RyueIdQtIINdu6XC7k5eVF3c/n86GxsZFVPNJNV5+8559/XlvSLMTlcmHBggWGBEXUVaF5qLjOLVHHKYqChoYGduzvRpIkwel0IhAIRLxeVVVUV1ezskodoivJ83q9mDt3Lg4dOgQA+OSTTzB37lycd955hgZH1BVmsxmnTp1q9QOFiNpWW1sb7xCSUqjZNpL6+nrIsszBFtQhun4SzJkzBzt27MD//u//oqioCC6XC4899hguuugio+Mj6hJJkrR1btnEQdQ+RVE4L16cCIIAv9+P+vp65ObmattDI3BZxaOO0r2sWV5eHsxmM2pqatCnTx8UFhYaGRdRTAiCAFEUceLECQSDwXiHQ5Tw6urqWC2KI7PZjLq6OsiyrG2rqamByaT765pIo+tV8+qrr2LFihW45557sHr1apx77rmYO3cuPv30U6PjI+oyQRAgCALXuSVqh6IoaGxsZBUvzkwmkza5u8fjQUtLC5M86hRdtd9Tp05h6dKlWvn4zjvvxGWXXYbVq1dj5MiRhgZIFAuhdW5PnjyJkpISViqIInA4HEwmEoDJZEJLSwvq6+s5GTV1ia538/z588P6BwDA4MGD8eyzzxoSFJERTCYTZFlGZWUl55kiOkswGITb7WaSlyAkSYLdbudnFXVJm+/mzZs3h13ev39/2OU33ngj9hERGUgURfh8Pq5zS3SWuro6JngJJiUlhU3n1CVtvqM3btwYdnn58uVhlz/66KPYR0RkMFEU4fF4OE0E0f8vGAyisbGRSR5RL9PmO7q9MjHLyNRTSZIEt9vd4YXBiXojh8PBihFRL9Rmktde53R2XqeeLDTDfENDQ7xDIYobWZbZF4+ol2pzyI6qqqitrdUqdpEux8q+ffuwfv16KIqCsrIyTJ06tVUs69evx969e2GxWDBz5kwMGDCgzWPfeOMNfPjhh9pC9bfddhsuvfTSmMVMPZ/ZbIbdbocoilznlpISq3hEvVebSZ7P58PDDz8ctu3sy7GgKAoqKirwy1/+EjabDfPnz0dpaSmKi4u1ffbu3Yvq6mqsXLkShw8fxrp167Bw4cJ2j50yZQquv/76mMdMvUdonduioiKkp6fHO5y4k2VZm5srNJm0yWQK+yeKojb/oMlk0v4O/aOeQZZlNDU1cTUYol6qzSRvw4YN3RLEkSNHUFhYiL59+wIARo0ahV27doUlebt378aYMWMgCAIGDRoEj8cDl8sFu93e7rFE7TGbzaiqqkJxcTEsFku8w+k2qqrC5/NpSZ3f70cwGIQgCJAkCaqqhv0LHXNmFf/spO7shC+UCJ59XWhbpATyzMtnH8NEMnZCVWwi6p0SYoZFp9MJm82mXbbZbDh8+HCrffLz88P2cTqd7R77/vvvY/v27RgwYADuuusuNslRVKIoapMlp6SkxDscQwSDQbS0tKCpqQk+nw+BQACqqrZKrEJilVC1tdLI2cnj2Ulk6O8z4wj9HamSGG372clkSkoKgsFg0iY5oYotq3hEvVdCJHmR+vad/cUSbZ+2jp00aRKmTZsG4HRV8tVXX8XMmTNb7b9161Zs3boVALB48eKwZNIIwWAQdru9U4mEKIpaH0P6j1idF1VV0djYiPPPP79XzDIfWuy8ubkZLS0tCAQCAE5XLntrIhtNKHkMBoPapNjBYBAWiwVpaWnIzc1Fenp60lQJjx8/jry8vIiPl58zkfG8RMdzE5koiobnFG1JiG8xm82Guro67XJdXR2sVmurfRwOR6t9ZFmOeuyZq3SUlZVhyZIlEe+/vLwc5eXl2uUz78cIiqLA7XZ36hd0dnY2GhsbDYiqZ4vleVFVFfv27UP//v171IhDRVHg8/m0Kl2o6TU3NxfNzc1h+7a0tMQpysSRnZ2tNVPX19fjxIkTEARBS/qys7N7bZXL7/ejqqoq6uPj50xkPC/R8dxEZrVaDc8pioqKol6XEN9gAwcORFVVFWprayHLMnbu3InS0tKwfUpLS7F9+3aoqopDhw4hPT0dVqu1zWNdLpd2/Oeff46SkpJufVzUM4UqxCdPnmyzmTHeQp3mq6ur8d133+Hbb7/FqVOn0NTUpDVDpqSk9IqKpNEEQYDZbIYkSQgGg2hoaMCxY8dw7NgxVFZWoqmpKaFfCx1lt9v5uiBKArrf5bIs4/Dhw3C5XBg1apRWCUhNTe1yEKIo4t5778WCBQugKArGjx+PkpISbNmyBcDpZtdhw4Zhz549mDVrFlJSUrRm12jHAsDrr7+OY8eOQRAEFBQUYMaMGV2OlZLDmevc9uvXL+5NeKqqIhAIoKmpSRsgIcsygNPz/YUGSlBsiKKo9dXz+/2orq4GcHqZqdTUVGRnZ8NiscT9ddEZPp8PXq+311Ypieg/BFXHZHffffcdlixZArPZjLq6Orz22mvYs2cPtm3bhjlz5nRHnN2qsrLS0NtXFAVHjx5lc20MGXVeZFlGeno6CgsLu/ULXVGUVgMkQtW5jg4U4Gsmss6eF0VREAwGYTKZYLFYkJGRgczMzB6TZJ86dQqBQKDN1zNfM5HxvETHcxOZ1Wo1vK9iW821uj6VXnrpJUyfPh1jxozBPffcAwAYPHgwXnjhhdhESJSgJElCc3Mz7HY7+vTpY9j9yLKM5uZmeDwerUqnqqpWpetMckfGOHMEsizLcDqdcDgckCQJqampyMrKQlpaWkL252QVjyi56EryTp48iauuuipsW2pqKvx+vyFBESWS0Dq3oiiGTdfTWW3NTReaZLinVIUovGk3VH0VBAEpKSnaAI5EGcnMvnhEyUXXu72goABHjx7FwIEDtW2hCYyJkoEkSaivr4coimGjtvVQFAVer1dL6tqbm456rtAADuA/o+hdLhckSQpr2o1HVdbn86GlpYVVPKIkoivJmz59OhYvXoyJEydClmX8+c9/xgcffID777/f6PiIEoYkSdo6n1lZWVH3CwQC8Hg8aG5u5gCJJBeadBk4/bpwOByw2+0wm83aAI7U1NRu6e9ZW1vL1x5RktH1jr/sssswf/58fPTRRxg8eDDsdjvmzp2LAQMGGB0fUUIxm82ora2FKIpIT0+HqqpoaWlp1fRqMpm0pI6VEwo5M8nyer1wu93dMjef1+uFz+fja5EoyehK8hobGzFgwAAmdUQ4/UUdmkg2EAhAURRIkqStu8oBEqTHmT8AQnPzOZ1ObTWS7OxspKenx6QpPzQwhIiSi653/cyZM/GDH/wAV155JYYPHx6TufGIejJJkrTRr0SxYNTcfM3NzaziESUpXd9Qa9aswaeffootW7bgpZdewqWXXoorr7wSw4YNY9WCiCjGzqzyqaoKj8eDxsbGTs3NxyoeUfLS9c7Pzs7G1VdfjauvvhoOhwMff/wx/vjHP2Lt2rWoqKgwOkYioqTW2bn5QoN/WMUjSk4d/nlXX1+P+vp6uN1uZGRkGBETERG1Qe/cfA6HgwkeURLTPRnyxx9/jE8++QR+vx8jR47EY489hvPPP9/o+IiIqA3R5uYTRZH9RomSnK53/xNPPIERI0ZgxowZGDJkSI9clJuIKBmcOTcfESU33WvX8tcgERERUc8RNXPbvn07xowZo/0dzYQJE2IfFRERERF1SdQk75NPPtGSvB07dkS9ASZ5RERERIknapI3f/587e8nn3yyW4IhIiIiotjQtV7OvHnzIm7/+c9/HtNgiIiIiCg2dCV5oeV1zqSqKmpqamIeEBERERF1XZtDZletWgXg9Azrob9D7HY7SkpKjIuMiIiIiDqtzSSvb9++Ef8WBAEXXnghRo4caVxkRERERNRpbSZ5N998MwDgggsuwCWXXNItARERERFR10VN8g4ePIjBgwef3kmScODAgYj7DRkyxJjIiIiIiKjToiZ5FRUVWLZsGQBg7dq1EfcRBKFVXz0iIiIiir+oSV4owQOA1atXd0swRERERBQbuqZQOduBAwfw9ddfxzoWIiIiIooRXUnek08+iX//+98AgE2bNuF3v/sdVqxYgbfeesvQ4IiIiIioc3QleSdOnMCgQYMAAB9++CGefPJJLFiwAB988IGhwRERERFR57Q5hUqIqqoA/rPyRXFxMQDA4/EYFBYRERERdYWuJO/CCy/Eyy+/DJfLheHDhwM4nfBlZWUZGhwRERERdY6u5toHH3wQ6enp+N73vodbbrkFAFBZWYnJkycbGhwRERERdY6uSl5WVhZuv/32sG2XXnqpIQERERERUdfpSvJkWcZbb72F7du3w+VywWq1YsyYMbjxxhshSbpugoiIiIi6ka4M7fXXX8c333yD++67DwUFBbDb7di4cSOam5tx9913GxwiEREREXWUriTvs88+w9KlS7WBFkVFRTjvvPPw2GOPMckjIiIiSkC6Bl6EplAhIiIiop5BVyVv5MiRWLJkCaZNm4b8/Hw4HA5s3LgRI0eONDo+IiIiIuoEXUneHXfcgY0bN6KiogIulwt5eXkYNWoUbrrpJqPjIyIiIqJO0JXkSZKE6dOnY/r06UbHQ0REREQx0GaSV1lZiTVr1uDEiRM477zzMHPmTPTp06e7YiMiIiKiTmpz4MXLL7+MPn36YPbs2cjLy8Mrr7zSTWERERERUVe0Wcn79ttvsXbtWqSkpGDw4MGYPXu2YYHs27cP69evh6IoKCsrw9SpU8OuV1UV69evx969e2GxWDBz5kwMGDCgzWObmpqwfPly2O12FBQUYM6cOcjMzDTsMRARERElijYrebIsIyUlBQCQmpoKv99vSBCKoqCiogKPP/44li9fjk8++QQnT54M22fv3r2orq7GypUrMWPGDKxbt67dYzdt2oShQ4di5cqVGDp0KDZt2mRI/ERERESJps1KXiAQwIYNG7TLfr8/7DKAmAzGOHLkCAoLC9G3b18AwKhRo7Br1y4UFxdr++zevRtjxoyBIAgYNGgQPB4PXC4X7HZ71GN37dqFp556CgAwduxYPPXUU7jjjju6HG9XqKoKr9cLr9eLQCCgzUGoqqr278zLiqKEXW5paYHb7Q7b/8xjzt4/0m2evX97MYQuh/bVu3+sCILQ7j7p6elobm6OyW3pFbotQRC0f2deNplMbV4f6XIsbvPsbX6/H01NTe3eT6TrO3IOI23vyPnuyO3G4v46ItprOtL2juxr5P21t19blyVJgsfjidntxeL4rn6uxOL8+3w+uN3umMbVW3i9XjQ1NXXomM68V2Lx/jLi/qMJBALIzs7u0DGx1GaSd+WVV6Kurk67PHr06LDLseJ0OmGz2bTLNpsNhw8fbrVPfn5+2D5Op7PNYxsaGmC1WgEAVqsVjY2NMY+9oxoaGvCDH/wg3mEQERGRwa666ir88Y9/jNv9t5nkzZw5s1uCiJQZn/2LPNo+eo5tz9atW7F161YAwOLFi8OSyVjLzMzEM888A7vdDkmSolZQQlWbs7eJoghVVSNWYiJVefRWbWK5b2hbLOj91SSKIoLBYEzuU4/2KqVt/Yu2r97t0aq10W5DEISwKqze2432uPVuN6py1dH7i7a9o58fsagoGlXZ1Ltve/cfut5kMkFRlHaP7+ztx+pyR3X1/IfOS6xvtzeI9n7qyu119LrOvO8683x05H6KiooMzSnao2uePKPZbLawCmFdXZ1WgTtzH4fD0WofWZajHpuTkwOXywWr1QqXyxW1ZFpeXo7y8nLt8pn3Y4S7774bR48ehdls7vCx2dnZCVGRTDQ8L9Hx3ETG8xIdz01kPC/R8dxEZrVaDc8pioqKol6na+1aow0cOBBVVVWora2FLMvYuXMnSktLw/YpLS3F9u3boaoqDh06hPT0dFit1jaPLS0txbZt2wAA27Ztw/Dhw7v9sRERERHFQ0JU8kRRxL333osFCxZAURSMHz8eJSUl2LJlCwBg0qRJGDZsGPbs2YNZs2YhJSVFa0qOdiwATJ06FcuXL8dHH32E/Px8PProo3F7jERERETdSVA5NKiVyspKQ29fURQ218YYz0t0PDeR8bxEx3MTGc9LdDw3kVmtVsNH17bVXBu1knfgwAFdNz5kyJCOR0REREREhoqa5K1duzbsstPphCAIyMrK0uZps9lsWLVqleFBEhEREVHHRE3yVq9erf391ltvoampCdOnT4fFYoHP58OGDRuQlZXVLUESERERUcfoGl37zjvv4Pbbb4fFYgEAWCwW3H777Xj77bcNDY6IiIiIOkdXkpeamoojR46Ebfvmm2+0pI+IiIiIEouuKVSmT5+OhQsX4rLLLtMmLt6zZw9++tOfGh0fEREREXWCriRvzJgxGDBgAD777DO4XC7069cPN910E4qLi42Oj4iIiIg6QfdkyMXFxbjxxhvR0NDQaskxIiIiIkosupI8j8eDdevW4bPPPoMkSXjttdewe/duHDlyBLfeeqvRMRIRERFRB+kaePHSSy8hPT0da9asgSSdzgsHDRqEnTt3GhocEREREXWOrkrel19+iRdeeEFL8IDTS5g0NDQYFhgRERERdZ6uSl56ejrcbnfYNofDwb55RERERAlKV5JXVlaGZcuW4cCBA1BVFYcOHcLq1asxceJEo+MjIiIiok7Q1Vx7ww03wGw2o6KiAsFgEGvXrkV5eTkmT55sdHxERERE1Am6kryGhgZMmTIFU6ZMCdteX1+P3NxcQwIjIiIios7T1Vw7e/bsiNvnzJkT02CIiIiIKDZ0JXmqqrba1tzcDJNJ1+FERERE1M3abK792c9+BgDw+/3a3yFNTU0YPXq0cZERERERUae1meQ9/PDDUFUVixYtwsMPPxx2XW5uLoqKigwNjoiIiIg6p80kb/DgwQCAiooKWCyWbgmIiIiIiLpOV6e65557Dl9//XXYtq+//hrLli0zJCgiIiIi6hpdSd7Bgwdx4YUXhm0bNGgQvvrqK0OCIiIiIqKu0ZXkmc1mtLS0hG1raWmBKIqGBEVEREREXaMrybv44ovx4osvorm5GcDp6VMqKipwySWXGBocEREREXWOrhUv7rrrLjz33HO45557kJWVhaamJlxyySWtRtwSERERUWLQleRlZmZi/vz5qK+vh8PhQH5+PpczIyIiIkpgupescLvd2L9/P7766ivk5ubC6XSirq7OyNiIiIiIqJN0j6595JFHsGPHDmzcuBEAUF1djZdeesnQ4IiIiIioc3Qlea+88goeeeQR/OIXv9BG1J5//vn45ptvDA2OiIg6TlEU+P1+BIPBeIdCRHGkK8mz2+0YOnRo2DZJkvgBQkSUAFRVRSAQQCAQgMlkQlZWFkpKSuIdFhHFma6BF8XFxdi3b1/YlClffvkl+vfvb1hgREQUXTAYhKIokCQJqampyMrKQlpaGkym//x2t1qtcDqdkCRdH/VE1MvoeuffeeedWLJkCYYNGwa/348XX3wR//rXv/DYY48ZHR8REeF0E2wwGITJZILFYkFOTg4yMzPbTOByc3NRX1/fjVESUSLRleQNGjQIS5cuxY4dO5Camor8/HwsXLgQNpvN6PiIiJKSqqqQZRmCIMBsNiMjIwPZ2dmwWCwQBEHXbQiCoFXzuEIRUfLRXcPPy8vD9ddfD7fbjaysLN0fMkS9kaqqAMD3AcXUmU2wFosF+fn5SE9PD2uC7aicnBy4XK4YRklEPYWuJM/j8eDll1/GZ599BlmWIUkSrrjiCtxzzz3IzMw0OkaihKKqKlRVRUZGBnw+H/x+PxRFgSiKrJZQh5xZrbNYLMjKykJWVhbMZnPM7kMQBOTl5cHhcLBvHlGS0fWOX7NmDUwmE5YsWYKCggLY7Xa88cYbWLNmDebNm2d0jEQJQ1VVKIqC/v37a1+YoZGNHo8HXq8Xfr8fsiwDOD0KndU+ClFVFcFgEKqqIiUlBWlpacjKykJqaqqhr5Ps7GxW84iSkK4k76uvvsKLL76IlJQUAKdH2z744IO4//77DQ2OKNEEg0EUFxeHVUQEQUBKSgpSUlJgtVq1/VpaWuDxeODz+RAIBKAoCkwmE6t9SUZRFMiyDLPZjJSUFGRlZSE9Pb1bXwehvnms5hElF13v9qKiItTW1qK4uFjb5nA4UFRUZFhgRIlGlmUUFRXBYrG0u68oisjIyEBGRgaA/1T7mpub0dzczGpfL3Z2E2xmZiays7Nj2gTbGazmESUfXUnekCFDsGDBAlx11VXIz8+Hw+HAjh07MGbMGHz00UfafhMmTDAsUKJ4CgQCOOecc5CWltap48+s9uXm5gI4XeHxer3weDzw+/3aCgUmk4nVlh5GlmWoqgqz2Yy0tDRkZ2cb3gTbUYIgwGazoaamJu4JJxF1D13fJIcPH0ZhYSEOHz6Mw4cPAwAKCwtx6NAhHDp0SNuPSR71RrIso0+fPlpVLlZMJlOrap8sy/B4PFq1L9R/SxTFLo2wpNgKzVkniiIsFgtyc3ORmZmZ8E3xmZmZqKuri3cYRNRNdCV5Tz75pGEBNDU1Yfny5bDb7SgoKMCcOXMijtjdt28f1q9fD0VRUFZWhqlTp7Z5fG1tLebMmaM1KV9wwQWYMWOGYY+DeqdAIACbzYbs7GzD7ys0H1pubm5Ytc/n86GpqSlsJK8gCBBFMaEqRb3ZmU2wKSkpyMzMRFZWltZPuacIVfNqa2tZLSZKArre5QcPHsTgwYNbbf/4449x5ZVXdimATZs2YejQoZg6dSo2bdqETZs24Y477gjbR1EUVFRU4Je//CVsNhvmz5+P0tJSFBcXt3l8YWEhli5d2qX4KHnJsozc3FxtMEU8mEwmpKWlhTUTh/r2eb1ebVAHAFb7YizUBNvWsmE9UVZWFpxOZ7zDIKJuoOvTatmyZXj99de1juIejwfLly/Hm2++2eUAdu3ahbFjxwIAxo4di127drXa58iRIygsLETfvn0hSRJGjRql7afneKKOCgaDyMzMRH5+frxDacVsNiMnJweFhYX43ve+hwEDBqBfv37IysqCKIpQFAWBQEBLUkgfRVG0JnKz2QybzYZzzz0X5557LgoLC5GRkdHjE7yQ/Px87fOciHovXZW8pUuXYu3atZg/fz6uvfZavPnmmxg2bBiWLFnS5QAaGhq0SonVakVjY2OrfZxOZ9gSajabTesb2NbxtbW1mDdvHtLS0nDrrbfi+9//fpfjpd4vGAwiNTUVffr0iXcoukSq9smyrI3kFQQBsixrKyn0lkSlq0IjngHAYrEgPT0d2dnZSElJ6fXN4BkZGWyuJUoCut7leXl5eOyxx/D444/jhRdewIQJEzrUv+2ZZ56JuEj2rbfequv4SNWI9j6ErVYr1qxZg6ysLBw9ehRLly7FsmXLkJ6e3mrfrVu3YuvWrQCAxYsXG169CQaDsNvtnerPI4pit/QP62lidV6CwSBSUlJw3nnn9ZovekmSEAgE4PV60djYqDXzBoNBCIKQ0FO4hFYXCf1TFCXs8+DMuEN/C4IAk8kEQRBa/X3mNrPZjHPOOadXVeg6IjU1FcePH4/4OcTPmch4XqLjuYlMFMW4tgjpSvKOHTuGlStXorCwELfddhteeeUVrFixAvfdd5+uEYdPPPFE1OtC6yparVa4XK6ILxKbzRY2Iqyurk6r3kU73mw2a9MEDBgwAH379kVVVRUGDhzY6vbLy8tRXl6uXXY4HO0+pq5QFAVut7tT0xhkZ2dHrHYmu1icl9BkxXl5eb1qBGJ+fr72eCRJ0pbOkmU5bAqXQCDQoZG8ZydgZ26LlDRGS8La+mcymbQpZUJ/h+I7e7/Q3x05Lw6HAz6fT/cxvY3P50NLS0ur7fyciYznJTqem8hCk5Abqa05i3Uleb/+9a9xxx13aFOk/OAHP8D69esxd+5crF27tkvBlZaWYtu2bZg6dSq2bduG4cOHt9pn4MCBqKqqQm1tLfLy8rBz507MmjWrzeMbGxuRmZkJk8mEmpoaVFVVoW/fvl2KlXqvUJJSXFycNFWdMxM+4PQ58Pl82iododGkQOvEDEBYwnX239GqZ5RY8vPzUVlZyXnziHopQdXRM7umpiZigrR7926UlpZ2KQC3243ly5fD4XAgPz8fjz76KDIzM+F0OvHCCy9g/vz5AIA9e/bg97//PRRFwfjx43HjjTe2efxnn32GN954Q/viufnmm3XHWllZ2aXH1B5FUXD06FFW8mKoK+cl0nq0vUmoYkXheF5OO3HihNZ0H8LPmch4XqLjuYnMarUa3ozdViWvzSSvvr5em68rkqNHj2LAgAFdiy4BMcnrebpyXmRZRklJSY+b80wvJjOR8byc5vV6cerUqbDPI37ORMbzEh3PTWTxTvLabD+ZPXt22OVQE2nI008/3YWwiOIvEAigX79+vTbBI2pPWloaLBYLp9sh6oXaTPLOftO73e42ryfqSWRZRlFREVJTU+MdClFcFRQUcN48gtU4SAAAH5JJREFUol6ozSSvvZFqiTrtAlF7AoEACgoKIk6pQ5RsUlNTkZqayh/uRL0Mh7tR0pFludvWoyXqKVjNSzxctYa6qs2hhD6fD08++aR2uaWlRbusqir8fr+x0RHFWCKsR0uUiCwWC9LS0rRVQCi+FEVBSkoKfD4fp7ihTmszyXvggQfCLo8fPz7scmjePKKeQJZlZGZmhi2RR0T/kZ+fjxMnTsQ7DMLpQkq/fv3gcDjQ1NQEURTjHRL1QG0meePGjeumMIiMJcsy0tPTe8x6tETxEKrmsYkwvgKBAPr27QuTyYT8/Hw0NTXFOyTqodgnj3q9YDAIi8WCwsJCDhYiakdBQQGbbONIURSkpaVpK9GYTCY+J9RpTPKoV1MUBZIkoaioiAkekQ4pKSnIzMxkNS9OFEVBYWFh2LasrCxWWKlTmORRr6WqKgRBSKr1aIliobCwkCNt40CWZeTn50fsf9e3b18+J9Rh/OajXim0Hi0TPKKOs1gsSE9PZ+WoG6mqCovFgpycnIjXS5IEq9WKYDDYzZFRT6ZrNXZZlvGPf/wDx44dQ0tLS9h1Dz30kCGBEXWWqqoIBoMoKSmBJOl6iRPRWfr06YNjx45x+o5uEgwGUVxc3OY+eXl5aGpq0lopiNqj6xtw1apVOH78OC677LKovzKIEkXow5Lr0RJ1niRJyMjIQEtLC6vhBpNlGXl5ee3+KBUEAX369MGpU6eYfJMuupK8L774AqtWrUJGRobR8RB1SSAQQFFRESwWS7xDIerxCgoKcOzYMSZ5BlJVFWazGbm5ubr2T0tLQ2ZmJrxeL58XapeuV0h+fj6Hb1PCC80txfVoiWJDkiRkZmZCUZR4h9JrybLc4emd+vTpw/6SpIuuSt6YMWOwdOlSXHvtta1+bQwZMsSQwIg6IhAIoKCgQJtbiohiIz8/H8ePH2fVyADBYBBWq7XDXUtCc+fV1tay3zG1Sder47333gMA/OEPfwjbLggCVq1aFfuoiDogNO0Av4SIYi/UN4/Ng7GlqipEUUReXl6njs/KykJDQwNkWeYgDIpKV5K3evVqo+Mg6hRZlpGVlYU+ffrA4XDEOxyiXqmgoADffvstk7wYCg0Q60qC1rdvX3z33Xes5lFUfMdSjyXLMjIyMrgeLZHBRFFEdnY2++bFiCzLyMnJ6fIAsdCADU6STNFETf/nzJmD5cuXAwB+9rOfRb2BtWvXxj4qonYEg0Gkpqaib9++8Q6FKCnYbDa43W5W87oo1Exrs9licnt5eXlwu90xuS3qfaImeffff7/298MPP9wtwRDpwfVoibqfKIrIysqCx+NhotcFsiyjX79+MfvsEgQBffv25dx5FFHUJO+iiy7S/h48eHC3BEPUHkVRYDKZutyXhYg6Lj8/n9W8LggGg8jOzkZaWlpMbzctLY0TV1NEfDVQj6GqKlRV5Xq0RHFiMpmQnZ3N9VM7SRAE5OfnG3Lbffv25dx51Aq/KalHOHM9WlEU4x0OUdKy2WxMJjohNFm7UT9QTSYTbDYbB2FQGCZ51CPIsozi4mL2OSGKM5PJhJycHFbzOkBRFGRmZhq+Gk9OTg7MZjOTcNK0m+QpioKnnnqKy5pR3AQCAfTr14/r0RIliM5O4JusVFXttqmeCgsLWc0jTbtJnslkQm1tLX8ZUFyE1nWMdUdlIuo8VvP0M7qZ9mwpKSmcO480ul5106ZNw0svvQS73Q5FUcL+ERklEAggPz8fmZmZ8Q6FiM5itVrjHULCUxQFGRkZyMjI6Nb7tdlsHJxGAHQua/bCCy8AALZv397qug0bNsQ2IiKcruDl5eUhJycn3qEQUQQmkwm5ublwuVxcVisKVVXjMmF7aO68yspK9mNOcrrematWrTI6DiKNLMvIzs5mvx+iBGe1WlFfXx/vMBJSdzfTni09PZ1z55G+JK+goADA6dJzQ0MDcnJy+KIhQwSDQWRkZGivOSJKXIIgwGq1wul0spp3BkVRkJaWhqysrLjG0adPHxw7dozf10lM17uyubkZL7/8Mj755BMoigJRFDFq1Cjce++9hg8Jp+QRDAZhsVi4Hi1RD5Kbm8tq3lkURUFhYWG8w9DWyK2rq2MSnqR0pffr169HS0sLli1bhtdffx3PPvss/H4/Xn75ZaPjoyShKArMZjPXoyXqYULVPI60PU2WZeTn5yfMpO25ubmcOy+J6Ury9u3bh4cffhhFRUXaF/HMmTPxxRdfGB0fJYHQerSxXLSbiLpPTk4O37s4PdDCYrEk3ICxwsJCJuFJSleSl5KSgsbGxrBtjY2NLP9Sl4V+XXI9WqKeK1TN+//au//YKO86DuDv58f9bHvX5577BbQd0LEpZotzJTjEsY6qiS6xaaZisulk/kgApxCmkoGpApMIyJwZjkVKyDQmYmJdYjInQdmE6SqIysgc3SbjR6E/ofxo4fnlH+SerOtd7wq9e5577v1Kluzunrt+nk+f+/Lp99dT6XuzGYbhimHa9/P7/bzncIUqqEq7//77sWHDBnzmM59BIpFAX18f/vCHP6ClpaXY8ZGHWZYF0zTR0NDgmqENIrox0Wi0oufmZbZ9cmvnRzwex6VLl5wOg0qsoKuxra0NiqLgwIEDGBwcRCwWw2c/+1k0NzcXOz7yMMMwUF9f79pGkYgKl+nN6+/vr7jvtGVZ8Pl8qK2tdTqUnARBQDKZRE9PD/fOqyB5v4mmaWLPnj1oa2vD/fffX4qYKoIkSdA0DZZlQRRFSJJUUXNaMvej9fv9TodCRFMkEolgaGjI6TBKTtd1NDQ0uL4Nr6qqQjgcxtWrVzk9pkIUdO/aP/7xjxxOm0KiKGLWrFmYPXs2brnlFiQSCYTDYfuvX13Xce3aNWia5slbx2mahmnTpvF+tEQeIwgCVFWFpmlOh1IyhmFAUZSy+YM1lUp58t8Vyq6gPvVFixbhT3/6Ez71qU8VO56KIooi/H6/PSk2w7Isu9AbGRnBtWvXYBgGNE2zi79M71+50XUdyWSy5PdyJKLSqK6uxsDAgNNhlIRlWZAkqazuzpOJlxtYV4aCfsPd3d148cUX8cILL0BV1TFd0j/4wQ9uKoBLly5h27Zt6OvrQyKRwMqVK7PekP7IkSPYtWsXTNPE4sWL0draCgB49dVXsWfPHpw+fRpPPvkkGhsb7ff87ne/w759+yCKIr7yla/gwx/+8E3FWiqCIMDn88Hn840rhmpra9HT02MXf5qmQdd1e9WUIAiuHfrVNA2qqo4paInIWzK9eb29vZ4vIgzDQF1dnSvb24nU1tbi4sWLME2z7GKnySnoG7h48WIsXry4KAF0dnbijjvuQGtrKzo7O9HZ2YmHHnpozDGmaWLnzp1Yu3YtVFXFmjVr0NTUhLq6OtTX12P16tV47rnnxrzn1KlTOHjwIH7yk59gaGgI69evx09/+tOyn4cgyzJCodC4oU7TNKFpGkZHRzE6OgpN02AYBnRdt7cpkWXZsS+0ruuora2FoiiO/HwiKp2amhoMDg46HUZR6bqOaDSKQCDgdCiTJggCUqkUTp48yUUYHlfQwotz586hra2tKBdDV1cX2tvbAVwfFm5vbx9X5HV3dyOdTtu3u1qwYAG6urpQV1eHurq6nJ+7YMEC+Hw+JJNJpNNpdHd347bbbpvyc3ADURQRCATGbcSZGfq9evUqRkZG7J4/XddhmiYsy4Isy0Utfg3DQHV1NeLxeNF+BhG5Szwex9mzZz3Zm5cZplVV1elQblggEEAkEsGlS5fKcuoPFSbvty+z8OJzn/tcUQK4cOGC3bujKMq4TZcBYHBwcMyXSVVVHD9+fMLPHRwcxJw5c+zHmTkIlea9Q7/vHwZ/77w/TdPGDP1mGjFRFG+q988wDASDQSSTyZs9FSIqI1VVVZ4s8IDrbacX7tATj8dx+fJlp8OgIirJwov169dn3SRzyZIlBb0/2z338n25JnOfvr1792Lv3r0AgE2bNrm6x0mW5aLHZxgGrl27hitXrmB0dHTM3D/LsiAIQkFDv4ZhwO/3Y9asWUVvDEuRl3LF3GTHvOQ2VbkJBAJ49913y2blaT6SJKGqqgqRSATTp093OpwpEQwGcerUqZseqZMkifOts5AkydF2piQLL9atW5fztWg0iqGhISiKgqGhoawXiaqqY1ZrDQwM5J3b9f73ZDZxzqalpWXM3Tv6+/sn/GwnxePxksaXGQIGrhfOmRW+V65csYd9M/v9AbB7/zL3o43FYiVZaVfqvJQT5iY75iW3qczN1atXMTo6OiWf5bRIJIKLFy8iFot56trRNA2XL1++qWk7kUgk60hcpctsEF5ME/3B4fjCi6amJuzfvx+tra3Yv38/5s2bN+6YxsZG9PT0oLe3F7FYDAcPHsRjjz2W93OffvppPPDAAxgaGkJPTw9uvfXWopxDpRAEwd7y5b1Dv5Zl2b1/mVW/pmli2rRpZb/QhYhuTjwex5kzZzwxwf/atWtIpVKea9dSqRROnDjhufOiAou8++67r2gBtLa2Ytu2bdi3bx/i8ThWrVoF4HrP244dO7BmzRpIkoSlS5di48aNME0Tzc3NqK+vBwC89tpr6OjowPDwMDZt2oSZM2fiiSeeQH19Pe655x6sWrUKoiji0Ucf5QVcJJnhW1mWEQ6HnQ6HiFwkHA7D7/eX/XYdpmmitra2LFfT5iPLMmKxGIaGhrgIw2MEa4LJax0dHVi6dKn9eN++fWNubbZlyxasXr26uBE64MyZM06HkBOHmLJjXnJjbrJjXnKb6tyMjIzg9OnTZd2bZ5om7r77bs8u4LMsC++++64973qyOFybnaIoRZ+rONFw7YRdW/v37x/z+Pnnnx/z+D//+c9NhEVERJUgFAohEAhMakGcm2iahmQy6enRIEEQkE6noeu606HQFJrwii3XLyQREblLIpEoywLCNE1UVVVVxK0YM3vnZe6gROVvwiKvnOdPEBGRewSDQQSDwbLrPLAsy96IvxJwWyFvmXDhhWEYOHr0qP3YNM1xj4mIiAqRSCTK6lZamqZ5cjXtRERRRCKRqIh7D1eCCX+D0WgUP//5z+3H1dXVYx5z40MiIipUIBBAMBiEruuuHykyTROhUAg1NTVOh1JyNTU1GB4ehqZprv890cQmLPKeeeaZUsVBREQVoFx680zTRDqddjoMx6RSKfzvf/9z/e+JJlY5fdBEROS4QCCAUCjk6rl5uq4jHo9X9J5xsixDURQuwihzLPKIiKik4vG4a1faWpaFQCCAaDTqdCiOi8ViEEXR1QU5TYxFHhERlZSbe/MMw6joYdr3EgQBqVTKtQU55ccij4iISs6N++bpuo5YLMZVpe8RCoVQXV3N3TTKFIs8IiIqOb/fj3A47JrePMuy4PP5UFtb63QorpNMJl3ze6LJYZFHRESOSCaTrunN03Ud6XSaW4ZkIYqiq+dRUm4s8oiIyBGyLKOqqsrxoUDDMKAoCvx+v6NxuFkkEinr+w9XKhZ5RETkmEQi4eg2HZZlQZIkxGIxx2IoF1yEUX5Y5BERkWNkWXZ0Yn9mNS2HafPz+XxQFIWFXhlhkUdERI6Kx+OO9Obpuo5oNIpAIFDyn12uMnvnUXngb4qIiBzlRG9eZphWVdWS/UwvEAQB6XQamqY5HQoVgEUeERE5rtRz83RdRyqV4jDtDQiFQq5YMEP5scgjIiLHSZKESCRSksLBMAxEIhGEQqGi/yyvSqVSXGlbBljkERGRK6iqWpIiL7PvG904URShqioXYbgcizwiInIFSZJQU1NT1EJP0zSkUikuHpgC0WgUPp+PPXouxquciIhcIx6PF63IM00T1dXVHKadQul0mr15LsYij4iIXEMURUQikaItwkgmk0X53Erl9/tRW1vLQs+lWOQREZGrqKo65UOAmqYhmUxymLYIVFWFJElOh0FZ8GonIiJXEUUR0Wh0ynrzTNNEVVUVqqqqpuTzaCxBEDBjxgzunedCLPKIiMh1pvJespZlIZVKTdnn0XiZIpp757kLizwiInKdqerN0zQNiUSCw7QlkEwmWeS5DK96IiJyJUVRbur9pmkiFAqhpqZmiiKiiWRuE8dFGO7BIo+IiFwp05t3o0WDaZpIp9NTHBVNpLa2lnvnuQiLPCIici1FUW7o/rK6riMej3PVpwPS6XRJ70NMubHIIyIi1xJFcdL7sFmWhUAggGg0WsTIKBe/31/UvQ6pcCzyiIjI1RRFmdTCCcMwOEzrMFVVb6gHlqYWizwiInI1QRCgKEpBvXm6riMWi0GW5RJERrmIoohkMsm98xzGIo+IiFwvGo3m7c2zLAs+nw+1tbUlioomUlVVhXA4zG1VHMQij4iIXC/TmzfRPC9d15FOpzlM6CKpVIpFnoNY5BERUVmYqDfPMAwoigK/31/iqGgikiQhFotx7zyHsMgjIqKykGtunmVZdjFB7lOJe+dZlgXDMBw/ZxZ5RERUNiKRyLi97zKraTlM606CICCVSnmyN880TWiaBk3T7POTJAmhUAixWAz19fWOxsflR0REVDYEQYCqqjh37hx8Ph90XUc0GkUgEHA6NJpAIBBAJBLBpUuXym6D6kyvnGmaEEURoihClmX4fD74fD6EQiH4/X5IkjTuDw2npw+wyCMiorJSXV2NgYEBe5hWVVWnQ6ICxONxXL582ekwsrIsC6Zp2otEMoVc5r9QKIRAIABZlie1Z6PTHC/yLl26hG3btqGvrw+JRAIrV65EdXX1uOOOHDmCXbt2wTRNLF68GK2trQCAV199FXv27MHp06fx5JNPorGxEQDQ29uLlStXYvr06QCAOXPm4Otf/3rpToyIiIoi05t35swZNDQ0cJi2TIiiiEQigbNnz8Ln8zkSg2ma9ly5TCEnSRJkWUYgEEAwGLR75bzA8SKvs7MTd9xxB1pbW9HZ2YnOzk489NBDY44xTRM7d+7E2rVroaoq1qxZg6amJtTV1aG+vh6rV6/Gc889N+6z0+k0Nm/eXKpTISKiEqmpqcGMGTMQCoWcDoUmobq6GqFQCJqmFa04f/+ih0wRJ8sy/H6/XcjJsuz5PxAc73Ps6urCokWLAACLFi1CV1fXuGO6u7uRTqeRSqUgyzIWLFhgH1dXV2f31hERUeXINupD7pdKpabkvraGYeDatWv2ogdBEODz+RAOh5FIJNDQ0IDZs2dj1qxZqK+vx7Rp06CqKqqqquDz+Txf4AEu6Mm7cOECFEUBcP3+hMPDw+OOGRwcHDPnQlVVHD9+PO9n9/b24jvf+Q5CoRCWLFmCD37wg1mP27t3L/bu3QsA2LRpE+Lx+I2cSknIsuzq+JzCvOTG3GTHvOTG3GTHvOQ22dzIsoy+vr68t5/LDK9mFj28d8FDMBgcM1fOjZy+ZkqSlfXr1+P8+fPjnl+yZElB78+2z0y+ClxRFGzfvh01NTV4++23sXnzZmzduhXhcHjcsS0tLWhpabEf9/f3FxSXE+LxuKvjcwrzkhtzkx3zkhtzkx3zkttkc2NZFq5cuWL/+54p5ARBgCiK8Pl8dkEXDoftQu79//aPjIxgZGRkSs9lKpXimploNLMkRd66detyvhaNRjE0NARFUTA0NIRIJDLuGFVVMTAwYD8eGBiwe/9yyVT6ADB79mykUin09PTYCzOIiIjIGYIgIJ1Oo6+vD36/31704PP5ymr1qts5nsmmpibs378fALB//37Mmzdv3DGNjY3o6elBb28vdF3HwYMH0dTUNOHnDg8P20uhz507h56eHqRSqak/ASIiIpq0QCCAuro6JJNJe69DFnhTy/FB7NbWVmzbtg379u1DPB7HqlWrAFyfh7djxw6sWbMGkiRh6dKl2LhxI0zTRHNzs72L9GuvvYaOjg4MDw9j06ZNmDlzJp544gkcO3YMv/nNbyBJEkRRxNe+9jVO0iUiIqKKIVhO31jNhc6cOeN0CDlxTkh2zEtuzE12zEtuzE12zEtuzE12Ts/JY78oERERkQexyCMiIiLyIBZ5RERERB7EIo+IiIjIg1jkEREREXkQizwiIiIiD2KRR0RERORBLPKIiIiIPIhFHhEREZEHscgjIiIi8iAWeUREREQexHvXEhEREXkQe/LKzPe+9z2nQ3Al5iU35iY75iU35iY75iU35iY7p/PCIo+IiIjIg1jkEREREXmQ1N7e3u50EDQ5s2fPdjoEV2JecmNusmNecmNusmNecmNusnMyL1x4QURERORBHK4lIiIi8iDZ6QDouiNHjmDXrl0wTROLFy9Ga2vrmNdfeeUV/P73vwcABINBfPWrX8XMmTMBAMuXL0cwGIQoipAkCZs2bSp1+EWVLzevv/46fvzjHyOZTAIA5s+fjwcffLCg95azfOf2wgsv4JVXXgEAmKaJU6dOYefOnaiurvb0NbN9+3YcPnwY0WgUW7duHfe6ZVnYtWsX/vnPfyIQCGDZsmX2cIqXrxcgf24qtZ3Jl5dKbWPy5aVS2xgA6O/vxzPPPIPz589DEAS0tLTg05/+9JhjXNHWWOQ4wzCsFStWWGfPnrU0TbNWr15tnTx5cswxb7zxhnXx4kXLsizr8OHD1po1a+zXli1bZl24cKGkMZdKIbk5evSo9aMf/eiG3luuJntuXV1dVnt7u/3Yy9fM66+/br311lvWqlWrsr5+6NAha+PGjZZpmtZ///tf+7vk5eslI19uKrWdyZeXSmxjLCt/Xt6rktoYy7KswcFB66233rIsy7KuXLliPfbYY+N+925oazhc6wLd3d1Ip9NIpVKQZRkLFixAV1fXmGNuv/12VFdXAwDmzJmDgYEBJ0ItuUJyU4z3ut1kz+3AgQP42Mc+VsIInTN37lz7u5LNP/7xD9x7770QBAG33XYbLl++jKGhIU9fLxn5clOp7Uy+vOTi9WtmMnmppDYGABRFsXvlQqEQZsyYgcHBwTHHuKGt4XCtCwwODkJVVfuxqqo4fvx4zuP37duHu+66a8xzGzduBAB84hOfQEtLS3ECdUChuXnzzTfx+OOPQ1EUPPzww6ivr590XsvJZM7t6tWrOHLkCB599NExz3v1mslncHAQ8XjcfqyqKgYHBz19vdyISmpnClFpbcxkVHob09vbi3feeQe33nrrmOfd0NawyHMBK8sCZ0EQsh579OhR/PnPf8YPf/hD+7n169cjFovhwoUL2LBhA6ZPn465c+cWLd5SKiQ3s2bNwvbt2xEMBnH48GFs3rwZTz/99KTyWm4mc26HDh0a00MDePuaySdX7rx8vUxWpbUz+VRiGzMZldzGjI6OYuvWrXjkkUcQDofHvOaGtobDtS6gquqYYZGBgQEoijLuuBMnTmDHjh14/PHHUVNTYz8fi8UAANFoFPPmzUN3d3fxgy6RQnITDocRDAYBAB/5yEdgGAaGh4cLzms5msy5HThwAAsXLhzznJevmXxUVUV/f7/9OJM7L18vk1GJ7Uw+ldjGTEaltjG6rmPr1q34+Mc/jvnz54973Q1tDYs8F2hsbERPTw96e3uh6zoOHjyIpqamMcf09/djy5YtWLFiBaZPn24/Pzo6ipGREfv///3vf6OhoaGk8RdTIbk5f/68/ZdRd3c3TNNETU1NQe8tV4We25UrV3Ds2LExr3n9msmnqakJL7/8MizLwptvvolwOAxFUTx9vRSqUtuZfCqxjSlUpbYxlmXh2WefxYwZM/DAAw9kPcYNbQ03Q3aJw4cPY/fu3TBNE83NzWhra8NLL70EAPjkJz+JZ599Fn//+9/t8f3MkvRz585hy5YtAADDMLBw4UK0tbU5dh7FkC83L774Il566SVIkgS/348vfelLuP3223O+1yvy5QUA/vKXv+DIkSP49re/bb/P69fMU089hWPHjuHixYuIRqP4/Oc/D13XAVzPi2VZ2LlzJ/71r3/B7/dj2bJlaGxsBODt6wXIn5tKbWfy5aVS25h8eQEqs40BgDfeeAPf//730dDQYA+1fvGLX7R77tzS1rDIIyIiIvIgDtcSEREReRCLPCIiIiIPYpFHRERE5EEs8oiIiIg8iEUeERERkQexyCMiukn9/f14+OGHYZqm06EQEdm4hQoR0Q1Yvnw5vvGNb+DOO+90OhQioqzYk0dERETkQezJIyKapJ/97Gf461//ClmWIYoiHnzwQfzqV7/Cr3/9a0iShPb2dnzgAx/A0aNHceLECXzoQx/C8uXLsWvXLhw6dAjTp0/HypUrkUwmAQCnT59GR0cH3n77bUQiEXzhC1/AggULHD5LIip37MkjIpqkb37zm4jH4/jud7+L559/Hvfcc8+4Yw4cOIAVK1Zgx44dOHfuHNauXYv77rsPHR0dmDFjBn77298CuH5vzw0bNmDhwoX4xS9+gW9961vYuXMnTp48WerTIiKPYZFHRFQEzc3NSKfTCIfDuOuuu5BKpXDnnXdCkiR89KMfxTvvvAPg+j0sE4kEmpubIUkSZs+ejfnz5+Nvf/ubw2dAROVOdjoAIiIvikaj9v/7/f5xj0dHRwEAfX19OH78OB555BH7dcMwcO+995YsViLyJhZ5REQOUlUVc+fOxbp165wOhYg8hsO1REQ3oLa2Fr29vTf9OXfffTd6enrw8ssvQ9d16LqO7u5unDp1agqiJKJKxp48IqIb0Nraio6ODvzyl79EW1vbDX9OKBTC2rVrsXv3buzevRuWZeGWW27Bl7/85SmMlogqEbdQISIiIvIgDtcSEREReRCLPCIiIiIPYpFHRERE5EEs8oiIiIg8iEUeERERkQexyCMiIiLyIBZ5RERERB7EIo+IiIjIg1jkEREREXnQ/wG4sO2sczdS1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Error plot in expected positive exposure\n",
    "up = np.mean(pi['tilde'], axis=1) -np.mean(pi['exact'], axis=1)+2.0*np.sqrt(np.mean(pi['tilde_var'], axis=1))   # 95% confident interval \n",
    "down = np.mean(pi['tilde'], axis=1) -np.mean(pi['exact'], axis=1)- 2.0*np.sqrt(np.mean(pi['tilde_var'], axis=1))   # 95% confident interval \n",
    "\n",
    "plt.figure(figsize = (10,6),facecolor='white', edgecolor='black')\n",
    "plt.plot(timegrid[1:], np.mean(pi['exact'], axis=1)-np.mean(pi['tilde'], axis=1), color = 'black', label = 'Error in Expected Positive Exposure')\n",
    "plt.fill_between(timegrid[1:].flatten(), np.array(down), np.array(up), color = 'grey', alpha=0.3)\n",
    "plt.grid(True)\n",
    "#plt.xlim(90, 110)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Error in Expected Positive Exposure')\n",
    "plt.legend(loc = 'best', prop={'size':10})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
